Nonlinear charge transport in DNA mediated 

by twist modes 



F Palmero a , JFR Archilla a ' x , D Hennig b and FR Romero 

a Nonlinear Physics Group (GFNL) of the University of Sevilla, Spain 

h Freie Universitat Berlin, 
Fachbereich Physik, Arnimallee 14, 1^195-Berlin, Germany 



Abstract 

Recent works on localized charge transport along DNA, based on a three-dimensional, 
tight-binding model (Eur. Phys. J. B 30:211, 2002; Phys. D 180:256, 2003), suggest 
that charge transport is mediated by the coupling of the radial and electron vari- 
ables. However, these works are based on a linear approximation of the distances 
among nucleotides, which forces for consistency the assumption that the parameter 
a, that describes the coupling between the transfer integral and the distance be- 
tween nucleotides is fairly small. In the present letter we perform an improvement 
of the model which allows larger values of a, showing that there exist two very 
different regimes. Particulary, for large a, the conclusions of the previous works 
are reversed and charge transport can be produced only by the coupling with the 
twisting modes. 
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1 Introduction 



Charge transport in DNA is an interesting subject because of its role for bio- 
logical functions such as DNA repair after radiation damage and biosynthesis, 
and, on the other hand, because of the possibility of building electronic de- 
vices based on biomaterials [1,2,3,4,5]. When a charge moves along the double 
DNA strand it is accompanied by a local deformation of the molecule, forming 
a polaron. Simple, but powerful models of DNA are based on the one proposed 
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by Peyrard-Bishop [6], in which the only variables are the distances between 
bases within each base pair. The corresponding hydrogen bonds are described 
by Morse potentials. In the context of polaron dynamics, a tight-binding sys- 
tem coupled with the bond distances is used for the description of the charge 
dynamics. In this the deformation of the helix is small a harmonic 

approximation of the bond potentials is considered enough giving rise to a 
Holstein model [7,8]. 

However, the helical structure of DNA is thought to play a key role in its 
functional processes and recently steric models of the molecule have been in- 
troduced in the context of vibrational motion [9,10,11,12,13] and charge trans- 
port [14,15,16]. As it will be shown in this article, the value of the parameter 
a that couples the electronic transfer integral and the distance between nu- 
cleotides is of particular importance, determining to a great extent the degree 
of charge localization and its mobility properties. However there is no exper- 
imental evidence of its value. Previous works [14,15,16] have dealt with very 
small values of a, because the angular deformations of the helix had to be 
small in compliance with the linear approximations performed in the expres- 
sions of the distances between nucleotides. One conclusion of these works are 
that charge transport in DNA is mediated by coupling of the charge carrying 
unit with the hydrogen-bond deformations, and that this transport survives 
to a small amount of parametric and structural disorder. However, we show in 
this work, taking into account the full equations of the model, that for larger 
values of a a completely different regime appears. In this regime charge trans- 
port is mediated by local winding and unwinding of the helix. Moreover, this 
regime is more robust with respect to disorder than the one mediated by the 
coupling with the hydrogen-bond distances. 



2 The model 



In order to describe the dynamic of a DNA chain, we consider a variant of the 
twist- opening model [9,10,11], which itself is a modification of the Peyrard- 
Bishop model taking into account the helical structure of the molecule and the 
torsional deformations induced by the opening of the base pairs. Each group 
of sugar, nucleotide and base will be considered as a single, non-deformable 
object, with mass m, which is an averaged estimation of the nucleotide mass. 
As we are interested in base pair vibrations and not in acoustic motions, 
we fix the center of mass of each base pair, i.e., the two bases in a base 
pair are constrained to move symmetrically with respect to the molecule axis. 
Moreover, the distances between two neighboring base pair planes shall be 
treated as fixed, because in the axial direction DNA seems less deformable 
than within the base pair planes [9]. 



2 




In our model, all bond potentials will be treated as harmonic ones. This can 
be justified because charge transport is related with small deformations of 
the double chain. Furthermore, as the angular twist and the radial vibrational 
motion evolve on two different time scales, they can be considered as decoupled 
degrees of freedom in the harmonic approximation [13]. Thus, the position of 
the n-th base pair is represented by the variables (r n , <f> n ), where r n represents 
the radial displacement of the base pair from the equilibrium value Rq, and 
<fi n its angular displacement from equilibrium angles with respect to a fixed 
external reference frame. A sketch of the model is shown in Fig. 1. 

The Hamiltonian of the system is given by H = H e i + H ra( i + H twist , where 
H e i corresponds to the part related with the charged particle described by a 
tight-binding system: H el = T,n E n\n){n\-V n - 1>n \n-l){n\-V n+ljn \n + l){n\, 
where \n) represents a localized state of the charge carrier at the n th base pair. 
The quantities {Ki.n-i} are the nearest-neighbor transfer integrals along base 
pairs, and {E n } are the energy on-site matrix elements. A general electronic 
state is given by \^/) = J2 n c n{t)\n), where c n (t) is the probability amplitude of 
finding the charged particle in the state \n). The time evolution of the {c n (t)} 
is obtained from the Schrodinger equation ih(d\^)/dt) = H e {\^) 

The nucleotides are large molecules and their molecular motions are slow com- 
pared to the one of a charged particle, the lattice oscillators may be treated 
classically and H ra d and H twist, describing the radial and torsional contribu- 
tions to H, are, de facto, classical Hamiltonians. They are given by (omitting 
the hat on them): 
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where p r n and are the conjugate momenta of the radial and angular coor- 
dinates, respectively. M = 2m and J = M it! 2 ,, are the mass and the inertia 
moment of each base pair respectively, Q r is the linear radial frequency, and 
is the linear twist frequency. We represent by Q n , n -\ — (<f>n ~ <f>n-i), the devia- 
tion of the relative angle between two adjacent base pairs from its equilibrium 
value 9 . 

The interaction between the electronic variables and the structure variables 
r n and <p n arises from the dependence of the matrix elements E n and 
on them. The first ones are given by [8] E n = E® + kr n , expressing the varia- 
tion of the on-site electronic energies {E®} with the radial deformations. We 
assume that the transfer matrix elements become smaller when the distances 
between nucleotides grows. That is, V njTl -i depends on the distances between 
two consecutive bases along a strand as V n ,n-\ = Vo(l — ad ntn -i), with d n) n-\ 
given by 



d n , n _i = [a 2 + (.Ro + r n f + (R + r n _i) 2 - 

2(^0 + r n )(Ro + r n _i) cos(^ + ^„,n-i)] 1/2 - k, (3) 

with l = (a 2 + 4_Rq sin 2 (^ /2)) 1//2 , and a is the vertical distance between 
base pairs. The parameter a describes the influence of the distances between 
nucleotides on the transfer integrals. In this paper we do not perform the linear 
approximation of the distances d n>n -i considered in Refs. [15,16], allowing for 
larger values of the perturbations. 

Realistic parameters for DNA are given in Refs. [10,17]. We have considered: 
a = 3.4A, m = 300 amu, #0 = 10A, fi r = 8x 10 12 s -1 , fi^ = 9x 10 11 s -1 , and 
Vo = 0.1 eV. We consider a and k as adjustable parameters. 

The present model is a modification of a previous one proposed in Refs. [14,15,16], 
in order to study the transport of charge by polarons in DNA. In these works, 
only changes on the relative angular velocity between two neighboring base 
pairs, ntn -i, are taken into account in the rotational kinetic energy, ignoring 
the contribution of the rest of the terms. Also, in the cited references, only 
very small displacements around the equilibrium position of the base pairs 
are considered, and the distances d n , n -i can be approximate by their Taylor 
expansion up to first order. In this paper we allow the spatial variables to have 
(relatively) large values and therefore we use the full nonlinear Eqs. (3) for 
the distances d„ „_i. 
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We scale the time according to t — > Q r t, and introduce the dimensionless 
quantities: f n = r n (M^/Uq) 1 / 2 , k n = k n / {M^Vofl 2 , E n = E n /V , Cl = 
Q^/Qr, V = V /(JQ 2 r ), a = a(V /Mn 2 y/ 2 , R = ^(M^/Vq) 1 / 2 . The 
scaled dynamical equations of the system, from which we have omitted the 
tildes, are: 



irc„ = (E n + kr n )c n 
-(1 - a d n+hn ) c n+ i - (1 - a d nn _i) c n _i , 



r n = 



a 
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(4) 
(5) 

(6) 



and the quantity r = hfl r /Vo determines the time scale separation between the 
fast electron motion and the slow bond vibrations. In the limit case of a = 
and uniform E n = Eq the set of coupled equations represents the Holstein 
system, widely used in studies of polaron dynamics in one-dimensional lattices, 
and for a = k = 0, and random E n , the Anderson model is obtained. Using 
the expectation value for the electronic contribution to the Hamiltonian, the 
new Hamiltonian H = (<f>\H\(f>) /Vq is given by 



h = E {^l + rl) + f w n + n 2 (<p n - ^f] + 

(£° + kr n )\c n \ 2 - (1 - arf n , n _i)(c*c„_i + CnC*^)^. (7) 

The values of the scaled parameters are r = 0.053, Q 2 = 0.013, V = 2.5 x 10~ 4 , 
Rq = 63.1 and Iq = 44.5. We fix the value k = 1 and consider the parameter a 
as adjustable, its value complying with the hypothesis of small deformations 
of the helix. 



3 Stationary polaron-like states 



As a first step in our study of nonlinear charge transport in this DNA model, 
we focus our interest in localized stationary solutions of Eqs. (4-6). Since 
the adiabaticity parameter r is small, the fastest variable are the {c n }, with 
a characteristic frequency (the linear frequency of the uncoupled system) of 
order 1/r ~ 19, followed by the {r n } with frequency unity, and the {</> n } 
with ~ 0.11. Using the Born-Oppenheimer approximation, we can suppose 



5 



initially that r n and <p n are constant in order to obtain the stationary localized 
solutions. For this purpose, we use a modification of the numerical method 
outlined in Refs. [8,18]. We substitute in Eq. (4) c n = $„exp(— iEt/r), with 
time-independent $„'s, and we obtain a nonlinear difference system E§ = 
A$, with $ = ($1, ...,<E»at), from which a map $' = is constructed, 

|.|| being the quadratic norm. 

Thus, using Eqs. (4-6), the stationary solutions must be attractors of the map: 



a 
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where d! = d(r',<j)'). In order to achieve convergence, it has been necessary 
to use the result of an interaction in the geometrical variables as a seed in 
the map corresponding to the probability amplitude. The starting point is a 
completely localized state given by c n = 5 nt o, r n = and <p n = 0, Vn. The map 
is applied until convergence is achieved. In this way both stationary solutions 
and their energies E are obtained. 

We analyze the ordered case, i.e., £° = E Q , that arises in synthetic DNA (the 
constant E can be set to zero by means of a gauge transformation), and the 
disordered case. In the latter, we introduce diagonal disorder in the on-site 
electronic energy E® by means of a random potential E® e [-AE, AE], with 
mean value zero and different interval sizes AE. Such deviations can be caused 
by inherent disorder in the surrounding water and by the inhomogeneous dis- 
tribution of counterions along the DNA duplex [19]. 

In the ordered case, we have found that in the ground state the charge is fairly 
localized at one site, and the amplitudes decay monotonically and exponen- 
tially with growing distance from the central site. The associated patterns of 
the static radial and relative angular displacements are similar. In the dis- 
ordered case, the localized excitation patterns do not change qualitatively. 
However, as the translational invariance is broken by the disorder, the lo- 
calized excitation pattern is not symmetric with respect to a lattice site, in 
opposition to the ordered case. 



In order to study the localization of the excitation, we introduce the degree 
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Fig. 2. Profiles of the stationary state for a = 0.05. (a) Ordered case, (b) Disordered 
case for AE = 0.01 eV. a) Electronic probabilities |c„| 2 . b) Radial displacements r n , 
c) Arc displacements Ro 9n,n-i- The displacements are in scaled units equivalent to 
~ 0.161. 

of localization L of an amplitude pattern {u n } as 



L = 



(EnM) 2 ' 



'ID 



and the participation number P as its inverse P = 1/L. This magnitude gives 
an estimate of the number of excited oscillators in the chain. 



We have studied the participation numbers of the electronic probabilities |$ ra | 2 , 
the radial displacements r n , the angular displacements 9 n>n -i, and the density 
of energy as functions of the parameter a. As Fig. 3 shows, for each given 
value of a and the disorder amplitude, the participation numbers of all these 
magnitudes are very close. When a increases, the participation decreases and 
therefore the localization increases. Comparing the ordered and the disordered 
case, we can observe that the localization is enhanced with the disorder, due 
to Anderson localization 



4 Charge transport in the absence of disorder 



In this section, we study charge transport along the double strand by moving 
polarons. Once a stationary state is obtained, it can be moved under cer- 
tain conditions. There exists a systematic way to do it, known as the pinning 
mode [21], which consists of perturbing the (zero) velocities of the ground 
state with localized, spatially antisymmetric modes obtained in the vicinity 
of a bifurcation. This method leads to moving entities with very low radi- 
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Fig. 3. Participation number P as a function of a for the electronic probabilities 
\$> n \ 2 (solid line), radial displacements r n (dotted line), angular displacements n ,n-i 
(dashed line) , and density of energy (dashdot line) of the ground state for the ordered 
case (AE = 0) and the disordered one (AE = 0.01 eV) 

ation but has the inconvenience of being applicable only in the neighbor- 
hood of certain values of the parameters. Instead, we use the discrete gradient 
method [22], perturbing the (zero) velocities of the stationary state {r n (0)}, 
{(/> n (0)} in a direction parallel to the vectors (Vr) n = (r n+1 — r n _i) and/or 
(V0) n = (0 n+ i — 4> n -i)- Although this method does not guarantee mobility, it 
nevertheless proves to be successful in a wide parameter range. 

4-1 Radial movability regime 

In absence of diagonal disorder, and if the parameter a is small enough 
(a < 0.005), it is not possible to move the polarons by perturbing the angu- 
lar variables. Mobility can only be accomplished through perturbation of the 
radial ones, as in Ref. [15]. This regime is explained in detail in this reference, 
but, for completeness, we summarize here the most relevant characteristics: 

(1) The charge probability moves with uniform velocity along the lattice with 
apparently constant profile. 

(2) The movement of the charge is accompanied by a relatively large radial 
deformation (r n ~ 0.02 A) and a negligible angular lattice deformation 

Wn,n-1~2-10- 5 A). 

(3) The velocity of the polaron increases approximately in a linear way with 
the kinetic energy of the perturbation, as shown in Fig. 6, until a critical 
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value, beyond which the localization is destroyed. 

If we observe the spatial variables, there exists a small radial amplitude os- 
cillation pinned at the starting site. There are radial and angular oscillations 
accompanying the charge and two larger angular oscillations, emerging from 
the same origin, that propagate in opposite directions with lower velocity than 
the polaron, as shown in Fig. 4. 

4-2 Twist movability regime 

However, for larger values of parameter a (a > 0.01), the polaron can only 
be moved by perturbing the angular variables. Perturbations of the radial 
variables cannot activate mobility. This is a consequence of the different con- 
tributions of the radial and angular components in the discrete gradient. 

On the twist movability regime we observe that: 

(1) The charge probability apparently moves with constant velocity and con- 
stant profile. 

(2) There are relatively large radial and angular deformations travelling with 
the charge ( r n ~ 0.051, R 9 n ,n-i ~ 2 A ). 

(3) The velocity of the polaron increases initially with the kinetic energy of 
the perturbation but tends to be approximately constant, as Fig. 6 shows. 
It decreases slowly with a. This moving excitation remains localized even 
for large perturbations. 

There remains a small radial oscillation localized at the starting site. Also 
an angular deformation travels with opposite direction and phase than the 
angular pulse moving with the charge. 

4-3 Mixed regime 

It is possible, for intermediates values of the parameter a (0.005 < a < 0.01), 
that the polaron can become mobile by perturbing any set of variables, the 
radial or the angular ones. Nevertheless, the movement is rather different as 
shown in Fig 6, when the polaron propagation is activated by means of radial 
perturbations, its characteristics are similar to the radial movability regimen, 
and likewise when it is activated by angular perturbations. 

In general, radial movability requires less energy and has higher velocity than 
the angular one, as shown in Fig. 6. There exists a small interval of pertur- 
bations energies, for which the movement can be activated by means of both 
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(a) 




Fig. 4. Radial movability regime. Polaron motion along the chain in ordered case 
for a = 0.0002 and an initial kick A r = 0.01 in the radial velocities, corresponding 
to a kinetic energy of the perturbation AK = 5 x 10~ 6 eV. (a) Uniformly moving 
electronic excitation, (b) Radial displacements (A), (c) Arc displacements i?o^n,n-i 
(A) with #0 = 10A. 
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(a) 




Fig. 5. Twist movability regime. Polaron motion along the chain in ordered case for 
a = 0.05 and an initial kick \ v = 0.01 in the angular velocities, corresponding to a 
kinetic energy of the perturbation AK = 0.02 eV. (a) Uniformly moving electronic 
excitation, (b) Radial displacements (A), (c) Arc displacements Ro9 njn -i(A) with 
R = lOA. 



11 



12 r 



10 



i 

i <P 
i 



03 



f 9 
i i 

<P 



^ — * — 



0.005 0.01 0.015 

A K (eV) 



0.02 



0.025 



Fig. 6. Velocity of the polaron as a function of the kinetic energy of the perturbation 
AK added to the ground state in the ordered case. Squares correspond to the 
radial movability regime (a = 0.0002). Circles correspond to the mixed regime with 
radial activation (a = 0.01). Triangles correspond to the mixed regime with angular 
activation for the same value of a = 0.01. Stars correspond to the twist movability 
regime (a = 0.05). 

kind of perturbations. The velocities obtained are similar but the evolution of 
the spatial variables are rather different. 

It is worth remarking that the limits of these regimes are not exact, depending 
on the kinetic energy of the perturbation, but we have intended to present a 
general picture. 



4-4 Tail analysis of the twist movability regime 



The consequences for the tail analysis of the radial movability regime has been 
analyzed in Ref . [16] . Here we use this simple but useful instrument to perform 
a similar study of the twist movability regime. The tails of the moving polaron 
have small amplitude and therefore their dynamics is to a great extent linear. 
The linear equations corresponding to Eqs. (4-6) in the ordered case are: 



irc„ = -c n+ i - c n _i ; r n = -r n ; n = -VL 2 (2 n - 0„_i - n+ i). (12) 

According to these equations, the variables r n are driven by the nonlinear 
terms and we will not consider them here. We propose the moving tail modes 
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c n = exp(— £/2 (n— vt) — i Et/r) and <f> n = exp(— £ (n—v t)), which will be valid 
(relatively) far from the head of the moving polaron and while (n — vt) > 0. 
These modes comply with the observed phenomena of similar exponential 
decay for |c n | 2 and <p n and that the second magnitude is a deformation and 
not an oscillation. Substitution into the first and third equations leads to: 

t£v = 4 sin(g) sinh(f/2), E = -2 cos(g) cosh(£/2), (13) 

v = 2Vt sinh^/2)^- 1 . (14) 

The velocity in physical units is v = a Q r v . The inverse characteristic length 
can be related to the participation number by 



(E"=_ooexp(-e 


M)) 2 


tanh(£) 


E"=_ 00 (exp(-e 


M)) 2 


tanh 2 (£/2) 



Certainly £ and P are determined by the full linear equations and depend on 
a, but in this way we can relate the values of £ with the numerical results 
for P. Within the twist movability regime the values of P change from 5 
to 12, which correspond to values for £, 0.74 and 0.33, and for v, 3.17 and 
3.11 (A/ps), respectively. For £ small, P ~ 4/£, v ~ and v = aQ r Q = 
3.4 A 8 ps _1 v / 0.013 ~ 3.1 A ps -1 . The numerical results vary from 2.1 to 3.4, 
i.e., they are in good agreement with the values of this tail analysis, taking 
into account the drastic simplification of the model. Moreover, they are also 
coherent with the numerical observation that the velocity of the polarons 
practically does not depend on the kinetic energy of the perturbation, as Fig. 6 
shows. The conclusion is that in this regime the polaron is basically riding on 
the angular wave. Equations (13) can be used to obtain the wave number q 
and the charge energy E and to compare them with the numerical data with 
good results, with q ~ 6 • 10 -6 and energies slightly below E = —2 (scaled 
units). This last value corresponds to the extended ground state with q = 
and £ = [16]. 



5 Charge transport in the presence of disorder 

If an amount of disorder in the on-site energies E® is introduced, with random 
values j-E^I < AE, we find that moving polarons exist below a critical value 
AE crit . Beyond this value polarons cannot be moved as found in Ref. [15]. The 
most relevant results are: 

(1) In the radial movability regime, mobility is very sensitive to disorder. 
Only a small degree of disorder allows activation of polaron motion (e.g., 
if a = 0.0002, AE crit w 0.005 eV). 
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(2) In the twist movability regime the moving polaron is very robust with 
respect to disorder (e.g., if a = 0.05, AE crit 0.05 eV). 

(3) For the mixed regime, angular activation leads to more robust (but slower) 
polarons than radial activation. If AE is high enough, only angular acti- 
vation is possible. The movement is very similar to the ordered case. 



6 Conclusions 

We have performed an analysis of the influence of the radial and angular 
perturbations on the properties of moving polarons in a three-dimensional, 
semi-classical, tight-binding model for DNA in both the ordered and disor- 
dered case. There exist three different regimes depending on the parameter a 
that describes the coupling of the transfer integrals with the deformations of 
the hydrogen bonds: a) Radial movability, with charge transport activated only 
by perturbing the radial variables; b) Twist movability, with charge transport 
activated by angular perturbations; c) Mixed regime. 

The properties of the moving polarons are different in each regime. In gen- 
eral, the mobility induced by angular activation is more robust with respect 
to parametric disorder, the polarons have lower velocity, and the activation 
energies are higher than in the radial movability regime. 

We have found that charge transport along the DNA chain is strongly in- 
fluenced by the value of a. Therefore, experimental observations of charge 
transport in real DNA could help to determine this crucial magnitude, which, 
in turn, is necessary to understand the phenomenon of charge transport in 
DNA. 



Acknowledgments 

The authors are grateful to partial support under the LOCNET EU network 
HPRN-CT-1999-00163. JFRA acknowledges DH and the Institut fur Theo- 
retische Physik for their warm hospitality 



References 

[1] M Ratner. Nature, 397(6719) :480, 1999. 

[2] HW Fink and C Schonenberger. Nature, 398(6726) :407, 1999. 



14 



[3] P Tran, B Alavi, and G Gruner. Phys. Rev. Lett, 85:1564, 2000. 

[4] E Braun, Y Eichen, U Sivan, and G Ben-Yoseph. Nature, 391(6669):775, 1998. 

[5] D Porath, A Bezryadin, S de Vries, and C Dekker. Nature, 403(6770) :635, 2000. 

[6] M Peyrard and AR Bishop. Phys. Rev. Lett, 62(23) :2755, 1989. 

[7] TD Holstein. Ann. Phys. NY, 8:325, 1959. 

[8] G Kalosakas, S Aubry, and GP Tsironis. Phys. Rev. B, 58(6):3094, 1998. 

[9] M Barbi. Localized Solutions in a Model of DNA Helicoidal Structure. PhD 
thesis, Universitat degli Studi di Firenze, 1998. 

[10] M Barbi, S Cocco, and M Peyrard. Phys. Lett A, 253(5-6) :358, 1999. 

[11] M Barbi, S Cocco, M Peyrard, and S Ruffo. Jou. Biol. Phys., 24:97, 1999. 

[12] S Cocco and R Monasson. Phys. Rev. Lett, 83(24):5178, 1999. 

[13] S Cocco and R Monasson. J. Chem. Phys., 112(22):10017, 2000. 

[14] D Hcnnig. Eur. Phys. J. B, 30:211, 2002. 

[15] D Hennig, JFR Archilla, and J Agarwal. Physica D, 180(3-4) :256, 2003. 

[16] JFR Archilla, D Hennig, and J Agarwal. In L Vazquez, MP Zorzano, and 
RS Mackay, editors, Proceedings of El Escorial 2002. World Scientific, 2003. 

[17] L Stryer. Biochemistry. Freeman, New York, 1995. 

[18] NK Voulgarakis and GP Tsironis. Phys. Rev. B, 63:14302, 2001. 

[19] RN Barnett, CL Cleveland, A Joy, U Landman, and GBSchuster. Science, 
294:567, 2001. 

[20] PW Anderson. Phys. Rev., 109:1492, 1958. 

[21] D Chen, S Aubry, and GP Tsironis. Phys. Rev. Lett, 112(23):4776, 1996. 
[22] M Ibahes, JM Sancho, and GP Tsironis. Phys. Rev. E, 65:041902, 2002. 



15 



